#ifndef PRIME_H_
#define PRIME_H_
#include <iostream>
#include <vector>
#include <map>
#include <string>
#include <gmp.h>
#include <gmpxx.h>
#include <random>
#include <time.h>
#include "math.hpp"

using namespace std;
/*
 * @brief 素数工具类
 */
class Prime{
	vector<int> primes; 	//小素数表
public:

	bool is_prime(int p){
		if(p == 1)	return false;
		if(p == 2) return true;
		int left = 2, right = sqrt(p);
		for(int i = left; i <= right; i++){
			if( p %i == 0) return false;
		}
		return true;
	}

	Prime(){
		//小素数表
		int tmp_primes[] = {
			2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,101,103,107,109,113,127,131,
			137,139,149,151,157,163,167,173,179,181,191,193,197,199,211,223,227,229,233,239,241,251,257,263,269,
			271,277,281,283,293,307,311,313,317,331,337,347,349,353,359,367,373,379,383,389,397,401,409,419,421,
			431,433,439,443,449,457,461,463,467,479,487,491,499,503,509,521,523,541,547,557,563,569,571,577,587,
			593,599,601,607,613,617,619,631,641,643,647,653,659,661,673,677,683,691,701,709,719,727,733,739,743,
			751,757,761,769,773,787,797,809,811,821,823,827,829,839,853,857,859,863,877,881,883,887,907,911,919,
			929,937,941,947,953,967,971,977,983,991,997,1009,1013,1019,1021,1031,1033,1039,1049,1051,1061,1063 
		};
		int pn = sizeof(tmp_primes) / sizeof(int);
		for(int i = 0; i < pn; i++){
			primes.push_back(tmp_primes[i]);
		}
	}
	
	vector<int>& get_primes(){
		return primes;
	}

	//素数检测，返回true表示n很有可能是一个素数
	bool miller_rabin(mpz_class n){
		if(n == 1) return false;
		if(n == 2) return true;
		mpz_class n_sub_one, d;
		n_sub_one = n - 1;
		d = n_sub_one;
		// n - 1 = 2^s * d
		int s = 0;
		while(d % 2 == 0){
			s++;
			d /= 2;
		} 
		default_random_engine e; 
		uniform_int_distribution<unsigned> u(0, primes.size() - 1);
		//随机筛选50个素数做检测
		int a;
		map<int, bool> a_mark;		
		mpz_class x, pre_x;
		for(int i = 0; i < 50; i++){ 
			for(int j = 0; j < 100; j++){
				a = primes[u(e)];
				if(a_mark.find(a) == a_mark.end()){
					a_mark[a] = true;
					break;
				}
			}
			x = Math::power_mod(a, d, n);
			pre_x = x;
			if(x == 1 || x == n_sub_one) continue;
			for(int r = 1; r <= s; r++){
				x = (x * x) % n;
				if(x == 1){
					if(pre_x != 1 && pre_x != n_sub_one){
						return false;
					}
					break;
				}
				pre_x = x;
			}
			if(x != 1){
				//一定为合数
				return false;
			}
		}
		return true;
	}

	
	/*
	 * 伪随机素数，返回0表示没有找到这样的素数
	 */
	mpz_class pseudo_random_prime(int bit_len){
		mpz_class a = Math::random_num(bit_len);
		default_random_engine e; 
		unsigned ts = (unsigned)time(0);
		e.seed(ts);
		uniform_int_distribution<unsigned> u(0, 9);
		u(e);
		string str;
		unsigned num;
		for(int i = 1; i <= bit_len; i++){
			if(i == 1){
				uniform_int_distribution<unsigned> ufirst(1, 9);
				num = ufirst(e);
			}else if(i == bit_len){
				num = 1;
			}else{
				num = u(e);
			}
			str += string(to_string(num));
		}
		//跳过5尾号为5的数字总是能被5整除
		int steps[] = {2, 4, 2, 2};
		mpz_class n(str);
		int k = 0;
		//限制搜索10万次
		int max_try_time = bit_len >= 1000 ? 100000 : bit_len * 100;
		for(int i = 0; i < max_try_time; i++){
			// cout << "n : " << n << endl;
			if(miller_rabin(n)) return n;
			n += steps[k];
			k = (k + 1) % 4;
		}
		return 0;
	}

};

#endif
